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Abstract. We use exact orbit integration to investigate particle accel¬ 
eration in a Gauss held proxy of magnetohydrodynamic (MHD) turbu¬ 
lence. Regions where the electric current exceeds a critical threshold are 
declared to be ‘dissipative’ and endowed with super-Dreicer electric held 
Efi = jyj. In this environment, test particles (electrons) are traced and 
their acceleration to relativistic energies is studied. As a main result 
we hud that acceleration mostly takes place within the dissipation re¬ 
gions, and that the momentum increments have heavy (non-Gaussian) 
tails, while the waiting times between the dissipation regions are ap¬ 
proximately exponentially distributed with intensity proportional to the 
particle velocity. No correlation between the momentum increment and 
the momentum itself is found. Our numerical results suggest an acceler¬ 
ation scenario with ballistic transport between independent ‘black box’ 
accelerators. 


1 Introduction 


Astrophysical high-energy particles manifest as cosmic rays or, indirectly, as 
radio waves. X-rays, Gamma rays. These often occur in transients, and with 
distinctly non-equilibrium energy distributions. A prominent source of sporadic 
radio- and X-ray emission is the Sun during the active phase of its 11-year cy¬ 
cle. Among the numerous mechanisms proposed for accelerating solar particles 
to high energies (see ^ for an overview), stochastic ones attracted particu¬ 
lar attention because they require generic input data and do not rely on spe¬ 
cial geometrical assumptions. In stochastic acceleration |2I3I4I5| . particles move 
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in random electromagnetic fields, where they become repeatedly deflected and, 
on average, accelerated. The electromagnetic fields are thought to arise from 
magneto-hydrodynamic (MHD) turbulence (e.g., 0 ), perhaps excited by the 
broadband echo of a magnetic collapse. The turbulence may host shocks and 
other forms of dissipation if critical velocities [7] or electric current densities 
uni are exceeded. Associated with dissipation are (collisional or anomalous) re¬ 
sistivity and non-conservative electric fields, which sustain, locally, the electric 
current against dissipative drag in order to meet the global constraints. However, 
a detailed balance on the level of individual charge carriers is impossible because 
the dissipative drag depends on particle position and velocity, whereas the elec¬ 
tric field is a function of position only. Thus the electric field may compensate 
the bulk drag, but a (high-energy) population can be left over and exposed to 
acceleration cm This lack of detailed balance is in the heart of dissipative accel¬ 
eration mechanisms. In plasmas, dissipation occurs at ‘ruptures’ of the magnetic 
structure, and is therefore localized around critical points of the magnetic field. 

The above scenario, first envisaged by Parker El for the solar atmosphere, 
has since been explored in a large number of numerical studies |12lldll4lbll5IH)ll7| 
On the theoretical side, most stochastic acceleration theories mm base on 
Fokker-Planck approaches, thus transferring two-point functions of the electro¬ 
magnetic fields into drift and diffusion coefficients of particles by probing the 
fields along unperturbed trajectories. Dynamical particle averages are then re¬ 
placed by field ensemble averages, neglecting the fact that real particles move 
in one realization of the random field. As a result, diffusive behaviour may be 
predicted even if particles are trapped in a single realization of the random field. 

In order to investigate the full diversity of orbit behaviour one must resort to 
numerical simulations. In the present contribution we analyze the behaviour of 
test particles in resistive MHD turbulence with localized dissipation regions, with 
particular emphasis on the validity of a Fokker-Planck description m- We use 
here exact orbit integration, and thus avoid any guiding centre approximations 
inE2- The price for rigorosity is computational cost, which makes the scheme 
only feasible with the aid of high-performance computing. 


2 Acceleration Environment 

The MHD turbulence has been been modeled by full 3D spectral MHD simula¬ 
tions and by Gauss field proxies |23I17| . We concentrate here on the latter, which 
is computed from the vector potential A(x, t) = ^j^.akCOs(k-x —a;(k)< —0k) by 
means of tabulated trigonometric calls. This allows to continuously determine 
the fields at the exact particle position, and avoids any real-space discretiza¬ 
tion artifacts, but the computational overhead restricts the k sum to a few 100 
Fourier modes ak. They are taken from the shell min(l“^)<|k|< with 

vl the rms thermal ion Larmor radius and li the outer scale of the power spectral 
density (|akp) oc {l + l'^k'^ + lyky + l1kl)~''. The electromagnetic fields are then 
obtained from 
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B = V X A 
E = -(9tA + 77(j)j 


( 1 ) 

( 2 ) 


where /ioj = V x B and r](j) = 77o6*(|j| — jc) is an anomalous resistivity 
switched on above the critical current jc ~ encg. Here, Cg and n are the sound 


speed and number density of the background plasma. The Gauss field A must 
satisfy the MHD constraints 


E-B = 0 if ? 7 (j) = 0 and E/B'^va 


( 3 ) 


with va the Alfven velocity. Equation can be achieved in several ways. 
For instance, one can use Euler potentials of which only one is time-dependent, 
or force A to point along a single direction. A somewhat more flexible way, used 
here, is axial gauge at • = 0 with dispersion relation aj(k) = k-v^. A constant 

magnetic field Bq along va can be included without violating Eq. ®, and we set 
= (Hq +ag)/{^onmp) with as the rms magnetic fluctuations and rrip the 
proton rest mass. In the present simulation, va and the background magnetic 
field are along the z direction. The total magnetic field B = y/B^ + a^ is a free 
parameter, which dehnes the scales of the particle orbits. In order to represent 
coronal turbulence we choose va ‘2 ■ 10® m/s, i' = 1.5, B ~ 10“^T, as = 
10“^T, Bq = 10“^T, and lx = ly = 10®m, G = lO'^m. The current threshold jc is 
exceeded in about 7% of the total volume. Note that our choice represents strong 
(cb Bq) and anisotropic (G ^ lx,ly) turbulence. To embed our simulation 
in the real solar atmosphere one should associate G with the radial direction in 
order to reproduce the predominant orientation of coronal filaments. 

3 Particle Dynamics 

3.1 Physical Scaling 

Time is measured in units of the (non-relativistic) gyro period f2~^=m/qB; ve¬ 
locity in units of the speed of light; distance in units of cf2~^. Particle momentum 
is measured in units of me; the vector potential in units of mc/q; the magnetic 
field in units of B; the electric current density in units of ilB/{^Qc), so that 
the dimensionless threshold current is // = {mjmp)cscfv\; and the electric field 
is measured in units of cB, so that the dimensionless Dreicer field is E'jj 
= {v'clT'){melm) with u/ the electron thermal velocity and t' the electron-ion 
collision time. The dimensionless equations of motion are 



dx' 


( 4 ) 


( 5 ) 
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Fig. 1. Evolution of electron kinetic momentum. Top panel: 200 sample orbits; 
adiabatic (a), and accelerated (b) cases. Bottom panel: electric current density 
along the orbits a) and b). The critical current density (|j| > jc) is marked by 
dotted line. The present simulation is an extension of the simulation of El 


with 7 the Lorentz factor, B' = V' x A', and j' = V' x B' the electric cur¬ 
rent. The dimensionless resistivity tj' is characterized by the resulting dissipative 
electric field Er? = 770 j relative to the Dreicer field Ed. We chose rj' such that 
Eq/Ed ~ 10 ^. 

3.2 Particle Initial Conditions 

We consider electrons as test particles. The initial positions are uniformly dis¬ 
tributed in space, and the velocities are from the tail u > 3 Vth of a maxwellian 
of 10® K, which is typical for the solar corona. Coulomb collisions are neglected, 
which is a good approximation once the acceleration has set on, but is not strictly 
correct in the beginning of the simulation. 

3.3 Numerical Implementation and Simulation Management 

Equations 0 and © are integrated by traditional leapfrog and Runge-Kutta 
schemes. The test particle code is written in FORTRAN 90/95 and compiled by 
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Fig. 2. Frequency distribution of the energy jumps Aj of Fig. ^ (black line), 
together with a best-fit Levy density (gray line) with parameters ao = 0.75, 
/3o = —0.26, and Cq = 0.035 (crosses). Inlets: cuts of the likelihood surface at 
(ao, Po, Co)- The 99% confidence level is marked boldface. 


the Portland Group’s Fortran 90 compiler (pgf90). Diagnostics and visualization 
uses IDL as a graphical back-end. The code is run on the MERLIN cluster of 
the Paul Scherrer Institut, and on the ANIC-2 cluster of the Universite libre 
de Bruxelles. The ANIC-2 cluster has 32 single Pentium IV nodes, a total of 
48 Gbyte memory, and Ethernet connections. The MERLIN cluster consist of 
56 mostly dual Athlon nodes with a total of 80 GByte memory, operated under 
Linux and connected by Myrinet and Ethernet links. MERLIN jobs are managed 
by the Load Sharing Facility (LSF) queueing system. Parallelization is done 
on a low level only, with different (and independent) test particles assigned 
to different CPU’s. MPICH/MPI is used to ensure crosstalk-free file I/O. The 
field data are computed on each CPU for the actual particle position. Random 
numbers are needed in the generation of the Fourier amplitudes and -phases of 
the electromagnetic fields, and in the particle initial data; they are taken from 
the intrinsic random number generator of pgf90. 
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Fig. 3. Left: travel times — tn between acceleration regions versus 

velocity Right: energy gain A'jn = 7n+i — In versus velocity. While 
scales with inverse velocity (solid line: best-fit), there is no clear trend in the 
energy gain Z\ 7 „. 


4 Diagnostics and Results 

In order to characterize the relativistic acceleration process we consider the evo¬ 
lution of the kinetic momentum P' = 7 v'. This quantity is directly incremented 
by the equation of motion 0 , and - ignoring quantum effects - can grow to ar¬ 
bitrarily large values, so that it can serve as a diagnostics of diffusive behaviour. 
Alternatively we may use the kinetic energy 7 . 

The results of the orbit simulations are shown in Figs. □-H When initially 
super-thermal {v ^ ivth) electrons move in the turbulent electromagnetic fields 
(Eqns.n] 0 , some of them may become stochastically accelerated. From a pop¬ 
ulation of 600 electrons we find that 35% of the particles are accelerated, while 
the other 65% remain adiabatic vmm during the simulation {fit < 8 ■ 10^). The 
two cases are illustrated in Fig. □ (top). The orbit a) conserves energy adiabati- 
cally during the whole simulation, while the orbit b) does not. The orbits of 200 
randomly chosen particles are also shown (gray) to trace out the full population. 
The bottom panel of Fig. □ shows the the electric current density along the or¬ 
bits a) and b). Time intervals where the critical current (doted line) is exceeded 
correspond to visits to the dissipation regions. As can be seen, acceleration (or 
deceleration) occurs predominantly within the dissipation regions. Accordingly, 
the orbit b) which never enters a dissipation region remains adiabatic. As a 
benchmark we have set rjo = 0 and found that no acceleration takes place at 
all, thus reproducing the ‘injection problem QH]’ Smaller 77 ' yield smaller (than 
35%) fractions of accelerated particles. 

A glance at graph a) of Fig. □shows that P'{t') is poorly represented by 
a Brownian motion m with continuous sample paths. Rather, P' changes 
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Fig. 4. Histogram of the quantity together with an exponential fit 

(dashed). 


intermittently and in large jumps. Indeed, if we consider the energy change 
^7n = 7n+i ~7n across a dissipation region, we find that its distribution P{A^n) 
has heavy tails and a convex shape which deviates from a Gaussian law (Fig. [3 
black line). In order to characterize P(Z\ 7 „) we have tried to fit it by a (skew) 
Levy stable distribution Pl{x). The latter is defined in terms of its Fourier 
transform m 


())l(s) =exp|-C'|s|“(l + z/3|^tan^)| (6) 

with 0 < a < 2, —1 < /3 < 1, and C > 0. The parameter a determines 
the asymptotic decay of Pl{x) ~ at x ^ (7^/“, and /3 determines its 

skewness. The probability density function (PDF) belonging to Eq. has the 
‘stability’ property that the sum of independent identically Levy distributed 
variates is Levy distributed as well. While rapidly converging series expansions 
of Pl{x) are available for 1 < a < 2 or large arguments x, the evalua¬ 
tion of Pl{x) at small arguments and a < 1 is more involved. We use here a 
strategy where Pl (x) is obtained from direct computation of the Fourier inverse 
Pl{x) = (27r)~^f e~^^^cj)L{s) ds, with the integrand split into regimes of differ¬ 
ent approximations. At small jsj, both the exponential in 4>l{s) (Eq.EJ and the 
Fourier factor are expanded; at larger |s|, 4>l{s) is piecewise expanded while 
g-isx jg retained. In both cases, the s-integration can be done analytically, and 
the pieces are summed numerically. The resulting (Poisson) maximum-likelihood 
estimates of the parameters (a,/3, C) are ao = 0.75, /3o = —0.26, Co = 0.035. 
The predicted frequencies are shown in Fig. [3 (gray line), and inlets represent 
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sections of constant likelihood in (a, /3, C')-space, with the 99% conhdence region 
enclosed by boldface line. The hnding ao < 2 agrees with the presence of large 
momentum jumps. 

In a next step we have investigated the waiting times Atn = tn+i—tn between 
subsequent encounters with the dissipation regions. Simple ballistic transport 
between randomly positioned dissipation regions would predict a PDF of the 
form f{\vn\Atn) with the particle velocity and f{x) the PDF of distances 
between (magnetically connected) dissipation regions. (There is no Jacobian 
d{At)/dx since we are dealing with discrete events.) This is in fact the case. 
FigurelSl(left) shows a scatter plot of the actual velocity Vn versus waiting time 
Atn- Gray crosses represent all simulated encounters with dissipation regions, 
including all particles and all simulated times. There is a clear trend for Atn 
to scale with Vn^, and the black solid line represents a best fit of the form 
Atn = L/vn with L = 9 • 20^ c/l7. When a similar scatter plot of velocity versus 
energy gain is created (Fig. fright), then no clear correlation is seen: the energy 
gain is apparently independent of energy. In this sense the dissipation regions 
‘erase’ the memory of the incoming particles. Returning to the waiting times, we 
may ask for the shape of the function f(\vn\Atn)- This can be determined from 
a histogram of the quantity \vn\Atn (Fig. 0 solid line). The decay is roughly 
exponential (dashed: best-fit), although the limited statistics does certainly not 
allow to exclude other forms. 


5 Summary and Discussion 

We have performed exact orbit integrations of electrons in a Gauss field proxy 
of MHD turbulence with super-Dreicer electric helds localized in dissipation re¬ 
gions. It was found that the electrons remain adiabatic (during the duration of 
the simulation) if no dissipation regions are encountered, and can become accel¬ 
erated if such are met. The resulting acceleration is intermittent and is not well 
described by a diffusion process, even if the underlying electromagnetic fields are 
Gaussian. On time scales which are large compared to the gyro time, the kinetic 
momentum performs a Levy flight rather than a classical Brownian motion. The 
net momentum increments in the dissipation regions are independent of the in¬ 
going momentum, and have heavy tails which may be approximated by a stable 
law of index 0.75. The waiting times between subsequent encounters with dissi¬ 
pation regions are approximately exponentially distributed, P{At) ^ ^-vAt/L^ 
indicating that the dissipation regions are randomly placed along the magnetic 
field lines. This one-dimensional Poisson behaviour is most likely caused by the 
Gauss field approximation, and the waiting times in true MHD turbulence are 
expected to behave differently. An ongoing study is devoted to these questions, 
and results will be reported elsewhere. In summary, our numerical results suggest 
that the acceleration process may be modeled by a continuous-time random walk 
with hnite or infinite mean waiting time, and infinite variance of the momentum 
increments. Such models can be described in terms of fractional versions [2T12 n) 
of the Fokker-Planck equation. 
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